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Abstract 

Morphological instability of the solid-liquid interface occuring in a crystal growing from an 
undercooled thin liquid being bounded on one side by a free surface and flowing down inclined 
plane is investigated by a linear stability analysis under shear flow. It is found that restoring forces 
due to gravity and surface tension is important factor for stabilization of the solid-liquid interface 
on long length scales. This is a new stabilizing effect different from the Gibbs-Thomson effect. A 
particular long wavelength mode of about 1 cm of wavy pattern observed on the surface of icicles 
covered with thin layer of flowing water is obtained from the dispersion relation including the effect 
of flow and restoring forces. 
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I. INTRODUCTION 



The interaction between an imposed shear flow and a phase transition underlies a broad 

rt 

range of phenomena lfl. The stability of order under the influence of shear flow is funda- 

q n 

mental for engineering to understanding frictional wear [2j and lubrication |3]. In pattern 
ormation in nature, ripple formation in sand induced by water shear flow are well known 
^J. In theoretical works, the effect of shear flow on the morphological stability has been 

studied flflQ 

An example of morphological instability of the solid-liquid interface in the long wavelength 
region of about 1 cm under shear flow bounded on one side by a free surface is wavy 
oattern occuring on the surface of icicles (see Fig. 1 in Ref. [8| and Fig. 9A in Ref. 

In its relevant experiment of a crystal growth from a thin liquid flowing down an 
inclined plane with angle 9 sketched in Fig. Q it is found that mean wavelength of the wavy 
pattern of the solid-liquid interface is about 0.83/(sin #)°- 6 ~ - 9 cm [10]. Ogawa and Furukawa 
have recently proposed a model of the mechanism of occurrence of the wavy pattern and 
obtained reasonable values of wavelength on the icicles However, in order to explain 
more quantitatively the experimental result mentioned above, we modify their dispersion 
relations in the form that includes the effect of restoring forces due to gravity and surface 
tension on stability of the solid-liquid interface. Furthermore, we improve their formulations 
by using a linear stability analysis under forced flow developed firstly by Delves [5|. From 
the dispersion relation in the long wavelength approximation, we present a new amplification 
rate and phase velocity for the fluctuation of the solid-liquid interface different from that of 
Ogawa and Furukawa's model. 

This paper is organized as follows. In Sec. II, we develop generally the dispersion rela- 
tion for the fluctuation of the solid-liquid interface. In Sec. Ill, a perturbed normal flow 
induced by deformation of the solid-liquid interface is determined in the long wavelength 
approximation. In Sec. IV, a general solution of the perturbed temperature distribution in 
the liquid is obtained. In Sec. V, we determine the dispersion relation for the fluctuation 
of the solid-liquid interface in a crystal growth from a thin liquid flowing down an inclined 
plane by applying the solutions in Sec. Ill and IV to the general formulation in Sec. II. 
Section. VI is devoted to discussion. Conclusion is given in Sec. VII. 
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FIG. 1: Schematic diagram of crystal growth from liquid flowing down inclined plane. 

II. DISPERSION RELATION FOR THE FLUCTUATION OF THE SOLID- 
LIQUID INTERFACE 

We consider a crystal growth from an undercooled thin liquid flowing down an inclined 
plane with the angle 9 10]. Hereafter the analysis will be restricted to two dimension in a 
vertical plane (x, y) sketched in Fig. ^ The primary shear flow is parallel to the x axis, and 
the mean velocity varies only in the y direction. The shear flow is bounded on one side by a 
free surface. We note that ho is the mean thickness of the liquid, and uq is the velocity at the 
free surface. In this section, we develop generally the dispersion relation for the fluctuation 



5, 



of the solid-liquid interface by following the ideas 

In the frame of reference moving at the solid-liquid interface velocity V, the equations 
for the temperature in the flowing liquid T} and that in the solid T s are 

&T X -dTi &T X &T X (d 2 T x d 2 T{\ 



dt dy dx dy \ dx 2 dy 



dT, T -3T S fd 2 T, d 2 T< 



dt dy s \ dx 2 dy 2 ) ' 
where t is time, u and v are the velocity components in the x and y direction measured in 
the laboratory frame in which the crystal is at rest, K\ and k s are the thermal diffusivities of 
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the liquid and solid, respectively. We substitute 7} = TJ + T/, T s = T s + T s ', u = U + u' and 

f = Ap\/ + v' into Eqs. (JTJ) and (0), where TJ, T s and U are unperturbed steady fields and 

T[, T s ', w' and v' are perturbed fields, respectively. Here ApV is the advection flow due to 

the density difference of the liquid and solid, Ap = (pi — p s )/pi, Pi and p s being the density 

of the liquid and solid. Then the equations for the unperturbed fields are 

d 2 % pVdfi n 

+ —-7^ = 0, (3) 



dy 2 k s dy 
and the equations for the perturbed fields are 



dy 2 Ki dy 



d 2 T, VdT s 

+ --7^=0, (4) 



dT! -dT' -dT! ,dft (d 2 T[ d 2 T'\ 



dt dy dx dy V dx 2 dy 



^ ^ =K .m + ^ } , (6) 



dt dy \ dx 2 dy 2 

where p = p s /p h 

Suppose that perturbations of the solid-liquid interface, temperature and normal flow are 
expressed in the following forms: 

((t, x) = ( k expfcrt + ikx], (7) 
T i = 9i(y) ex P exp [at + ikx), (8) 

T 's = 9s{y) exp (-T^VJ exp[(xt + ikx], (9) 

v' — Vk expfai + ikx], (10) 

where k is the wave number, and a = a r +i&i, u T being the rate of amplification or damping 
and —Oijk being the phase velocity of the disturbance, (k and Vk are the amplitudes of 
perturbed interface and perturbed normal flow, gi and g s are the amplitudes of perturbed 
temperature of the liquid and solid, respectively. Substitutions of them into Eqs. © and 
© yield 

^ iu2 ■ (P V \ 2 ■ ° ■ MU(y) \ v h d% fpV \ 
\ k +[■*-) + — + } 9i = — -77 exp — y , (11) 



dy 2 | \2kJ «j «/ J ' «j V 2ft; « 
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d2 9s J , 2 . (PV 
dy 2 r \2k. 



+ 



a 



9s = 0. 



(12) 



The following calculations are to first order only in the amplitude of the initial perturba- 
tion. The continuity of the temperature at the perturbed solid-liquid interface, y = ((t,x), 



is 



(T, + T()\ y=c = (f s + T%=C = T m + G(k)C, (13) 

where T m is the melting temperature, and G(k)( is the temperature difference from T m 
due to a deformation of the solid-liquid interface. The form of G(k) will be specified later. 
Linearizing Eq. (fT^j) at y — 0, Eq. (fT^j) gives to the zeroth order in (k, 



Ti 



i\y=o 



s\y=0 



T 

J. rr 



and to the first order in (k, 



dy 



(k + gi 



y=0 



dTs 
dy 



y=o' 



Ck + g 



s\y=0 



G(k)Ck- 



It follows from Eq. (fTK|) that the amplitudes of gi\ y =o and g s \ y =o are of order (k- 



gi\ y =o 



dJ\ 
dy 



+ G(k) C 



(14) 



(15) 



(16) 



gs\y=0 



dT a 



dy 



y=0 



+ G(k) C 



The heat conservation at the perturbed solid-liquid interface is 



d(f s + T' s ) 



dy 



K, 



dm + T{) 



y=< 



dy 



y=< 



(17) 



(18) 



where L is the latent heat per unit volume and K s and Ki are the thermal conductivities of 
the solid and liquid, respectively. We linearize in the same way Eq. (J18|) at y — 0, Eq. (|18|) 
gives to the zeroth order in 



LV = K, 



dTs 

dy 



y =o dy 



y=o 



(19) 



and to the first order in 

La( k = K, 



d 2 T s 



dy 2 
d 2 Ti 



dy 2 



t P v i , dg s 

y=o 2k s dy 

pV dgi 

y=o 2ki dy 



y=0 



y=o 



(20) 



By substitutions of Eqs. (fTfij) and (|T7|) into Eq. (j20J), the dispersion relation for the fluctu- 
ation of the solid-liquid interface becomes 



(d 2 T« 



a 



L 



dy 2 



y=0 



Kj, 
L 



d 2 T, 



dy 2 



+ 



EL. 
2ki 



- Q, 



dT« 



dy 
_dJ\ 
dy 



y=0 



y=o 



G(k) 



+ G(k) 



(21) 



where we have defined the so-called propagator in the liquid and solid as follows 



Qi 



dsn 

dy 



y=o 



Q, 



dgs_ 
dy 



y=0 



(22) 



9l\y=0 9s\y=0 

which describe the motion of the interface in response to the propagation of a temperature 
disturbance, here it is the latent heat release. The general formulation above will be applied 
in Sec. V. 



III. THE PERTURBED NORMAL FLOW OVER THE SINUSOIDAL INTER- 
FACE 

It is firstly necessary to know the primary shear flow field U(y) and the amplitude Vk 
of the perturbed normal flow in Eq. (|11)1. We determine the perturbed normal flow over 
the interface to first order only in the amplitude of initial perturbation by following the 
formulation of Benjamin Q|. In his treatment, the bottom is flat interface and does not 
move. When the crystal grows, however, the solid-liquid interface moves and it may be 
not flat if a morphological instability occurs. Therefore, it is necessary to modify several 
boundary conditions in his formulation. We make use of the same dimensionless variables as 
those used in Benjamin's paper, which are different from those used in Ogawa-Furukawa's 
paper ^]. Hereafter we refer to their models as O-F model. 

With reference to Fig. Q the primary shear flow assumed steady, is pararell to the x axis, 
with the velocity U varing only with y. In the frame of reference moving at the solid-liquid 
interface velocity V, the Navier-Stokes equations are 

du ydu du du 1 dp ( d 2 u d 2 u\ . ^ 
dt dy U dx V dy p t dx \ dx 2 dy 2 y ^ s ' 

dv - dv dv dv 1 dp ( d 2 v d 2 v\ 
dt dy dx dy p\dy \dx l dy 2 J 
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where t is time, u (v) the velocity components in the x (y) direction, p the pressure, pi the 
liquid density, g the gravitational acceleration, v the kinematic viscosity, and 9 the angle of 
the inclined plane. The equation of continuity is 

du dv q 
dx dy 

In this section the coordinates (x, y) and velocities (u, v) are made non-dimensional by 
taking the mean thickness h of flowing liquid as the unit of length and the velocity uq at 
the free surface as the unit of velocity respectively. By substitutions of = (x,y)/h , 

(u*,v*) = (u,v)/uq, p* = p/(piUq) and t* = tuo/ho, the equations of motion and continuity 
can be written in the following dimensionless forms: 

du* -du* du* du* dp* 1 ( d 2 u* d 2 u*\ sinO 

- V *7^7+ u *7^r + v *7^7 = -7^T + Tr:[l^ + ^)+^r' ( 26 ) 



dt* dy* dx* dy* dx* Re V dx 2 dy 2 J F 2 



dv* T - r dv* dv* dv* dp* 1 ( d 2 v * d 2 v*\ cos 9 



dt* dy* dx* dy* dy* Re \ dx\ dy 2 J F 2 

^ + ^ = 0, (28) 
dx* dy* 

where Re = uoh^/v is the Reynolds number and F = Uq/ '(gho) 1 ^ 2 is the Froude number. 
Let 

u* = U* + u'*, v* = ApV* + v'*, p* = P* + p'*, (29) 

where U* and P* are the dimensionless velocity and pressure of the primary flow, and primed 
quantities are the dimensionless velocity and pressure perturbations. Substitution of Eq. 
flU into Eqs. lj2E j) . p7 |l and yields 

1 d 2 U* - dU* 1 

* + P V*-^ + — sm6 = 0, (30) 



Re dy 2 dy* V- 

dP* 1 



dy F 2 



cos# = 0, (31) 



du'* - du'* - d< dC/» , 9pi 1 / d 2 u'* d 2 u'* , , 



<9t* <9y* <9x* dy* * cte* Re \ cte 2 ch/ : 

9v* y® v * fj ® v * ®P* 1 f9 2 v* d 2 v^ . , ^ 

<9£* * dy* * dx* dy* Re V <9x 2 <9y 2 ' 
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^ + ^ = 0, (34) 
or* ay* 

if quadratic terms in the perturbation quantities are neglected. Using the typical values of 
V ~ 10~ 6 m/s and «o ~ 10™ 2 m/s in the experiments P, Oh we can neglect the pV* term 
in Eqs. (|3T)j). (jH^j) and because V* is the ratio of V to «o- 
Under the boundary conditions, 

U* = o (y* = o), ^ = (y* = l), P* = Po* (y. = i), (35) 

ay* 

the solutions of Eqs. (J37IJ) and (f3~Tj) are respectively jisl ]. 

& = 2y* - y 2 , (36) 



P* = P * + ^(i-y*), (37) 

where Po* is the dimensionless pressure of atmosphere. Eq. (jMj) allows the use of a stream 
function if)', in terms of which u'^ and v'^ can be expressed as follows: 

/ , dip' 

ay* or* 

Eqs. (jH2I) and can then be written as 

dt^dy* * dx^dy* dy* dx* dx* Re \dx1dy* dyf J ' 

9V - <9V ap' 1 /a 3 V>' aV \ , < 

+ U*—^r = ^ + —[—^ + - . (40) 



dt*dx* dxl <9y* Re \ ax 3 dx*dyl 
If the perturbation of the solid-liquid interface is represented in a dimensionless form, 

C*(**,£*) = <5bexp[cx*t* -M«x*], (41) 

the corresponding perturbations of the stream function, pressure and the liquid-air surface 
may be written as respectively, 

i 1 ' = $bf(y*) exp[cr*£* + ifix*], (42) 



p[ = 5bU(y*) exp[cr*t* + ipix*), (43) 
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= 1 + 5 t exp[cr*i* +ifxx*), (44) 

in which Sb = Cfc/^o an d St = ik/ho are dimensionless amplitudes of the solid-liquid interface 
and the liquid-air surface respectively, \i = kho is the dimensionless wave number, and 
cr* = crho/uo- When we substitute Eqs. and (jUJJ) into Eqs. (JUS) and (jUJ, and II is 

eliminated from them by cross differentiation, the linearized equations of motion lead to the 
Orr-Sommerfeld equation: 

d 4 f 9 2 d 2 f , 4f 



dyi dyl 



The perturbed flow was assumed to be stationary from the outset in the O-F model. This 
formally amounts to neglecting the cr*/ fi term in Eq. (|45)1 . This assumption will be justified 
in Sec. VI. Since the value of the mean thickness ho is about 10~ 4 m, and the typical value 
of wavelength of the wavy pattern observed on the surface of icicles is about 1 cm ^, Q| > 
the value of /i = kho is about 6 x 10~ 2 |8j. Therefore, in the long wavelength approximation 
retaining up to the first order in /i, Eq. (|45jl becomes 



f// ¥ Re{(2y*-yt)Pi + 2f\, (46) 



dyf [ * dyl 

where we have substituted Eq. for U m . We note that Re becomes 0(1) when we use 
the typical values of uq and ho used above, and v = 1.8 X 10~ 6 m 2 /s of water, therefore the 
primary shear flow is laminar j^, 0, . 

The problem entails five boundary conditions as follows. Since both velocity components 
must vanish at the perturbed solid-liquid interface, we have 

<L=C*-Ap|i = 0, (47) 



dC 

(£/* + <)U =f , + ApV*^ = 0. (48) 
The kinematic condition at the free surface is 

^L + fj f?k = v ' I , f49) 
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At the free surface the shear stress must vanish and the normal stress must just balance the 
normal stress induced by surface tension; 



du* 



+ 



y*=£* dx* 



0. 



-p*\ v .=e, + 



2 dv* 



Re dy* 



O r, — U, 



(50) 



(51) 



where S=7/ (pih ul), 7 being the surface tension of the liquid-air. Linearizing Eqs. (jUj) and 
d3Hj) at = and Eqs. {EJ), and (JS3J at = 1, Eqs. (JEJ) to (jHTJ) become respectively, 



(7* 



f\ y ,=o = iAp— , 
p 



dy* 



y*=o 



-2 — ipApV*. 



. ///Re cos o_, ,., , . 
5 b - 1 [ ^—^ h /TReS S t 



F 2 



ipRe I J +3/i 2 l ^ 



d 2 / 

5 6 - pV*Re—j 

w.=i 



4. 



If we formaly put 



/(y«) = 



AT 



N=0 



(52) 
(53) 
(54) 
(55) 



(56) 



(57) 



then this series is seen to constitute a solution of Eq. (|4l)Jl when the coefficients An are 
made to satisfy the following recursion relation 

N(N - 1)(N - 2)(N - 3)A N 

= 2ipRe(N -3)(N- 4)^_ 3 + ipRe{2 - (N - 4)(iV - 5)}A N ^, (58) 

for iV > 3. Eq. (|58|) gives every other A^ in terms of the first four coefficients A to A3. 
The approximation to the series solution up to the first order in p requires seven coefficients 
of the expansion Eq. (|57|). Therefore, the other coefficients are given as follows: 



^4 = 0, 



(59) 
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(60) 



A e = (61) 



A 7 = -^A 3 . (62) 
7 210 3 y J 



Hence, the approximate series solution can be written as 

ifjRe 4 \ / i/xRe 5 



/(».) = ( 1 + -^j-K I A, + ( !/• + "^-K ) 



+ (y> + ^) A 2 + (rf + J^rf - |^Jj (63) 

The four constants A to A 3 of the solution of the fourth order Eq. ()46|) are determined 
from the boundary conditions Eqs. f|52(l to (f56j) in the form neglecting the terms including 
K, cr^/fj, and p 2 . 

First, the boundary conditions Eqs. (JH2J and (JHHJ) give respectively, 

= 0, (64) 



A l = -2. (65) 
Eliminations of S t from Eqs. (|54*jl and and from Eqs. (J54*j) and (J5fi|) yield respectively, 



dyl 



-2/| 



(66) 



ci 3 / 



where 



dyl 

fiRe cos 6* 



2/iHe- — 

dy* 



!/*=! 



n g fin COS 9 1 ^hl , 
+ /TReS = — fc + -i-5i-fc . 



(67) 

(68) 
. When we 



F 2 I/M p^z/Mo 

Eq . O represents the restoring forces due to gravity and surface tension IllbL 
use the typical values of uq ~ 10 2 m/s and /io ~ 10 4 m in the experiments |2llld], and the 
physical properties of water, p\ = 1.0 x 10 3 Kg/m 3 , v — 1.8 x 10~ 6 m 2 /s and 7 = 7.6 x 10 -2 
N/m, a becomes 0(1) for the wavelength of wavy pattern occuring on the icicles and the 
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inclined plane, therefore we treat a as zeroth order in terms of fi in the following calculations. 
Since uq and h Q are not independent quantities, this rough order estimate will be justified 
more quantitatively in Sec. VI by using other parameters being controlled easily in the 
actual experiment. 

Using Eqs. (pjH|) . (jMj) and (p3o"|) . Eqs. and (Jo"7)) give respectively, 

4+ ll^e\ 2+ f 8 + M?te\ A 3=4 + 1^, (69) 



-(l + ^V 2 +f- + 6-^V 3 = -f2 + ^). (70) 



30 J V 35 J ° V 30 

Retaining up to the first order in /i, the solutions of these simultaneous equations for A 2 
and A 3 are expressed as follows: 

3(2 — ia) „ —96 — 8ia , , 

6 — «a 105(6 — 2a) 2 

^ 4za , , 

^ = ^+" Rea 35(63l5F' (72) 

When these expressions of A to A 3 are substituted into Eq. (jHSJ), the final form up to the 
first order in /x is 



3(2 — ia) 2 ia , f —96 — 8ia 2 4za , 

f{y*) = -2y* + p — : — V* + 7, — -y* + ^a < — — — -y, + — — -y, 

6 — ia 6 — ia y 105(6 — lay 35(6 — iaY 



+ 15(6 -laf* 30(6 - iaf* + 210(6 - taf*} ' (73) 

Applying this result at = 1 to Eq. (jSIJ), we can know the relation between the amplitude 
and the phase of perturbation of the solid-liquid interface and that of the liquid-air surface: 



St = -f\»=i6 b . (74) 
In the article 8j, the following function was obtained: 

f(y*) = -2y* + yl (75) 

If we substitute Eq. (175)) into Eq. (JHJ), 8 t = <5&, which indicates that the liquid-air surface 
fluctuates with the same amplitude as the solid-liquid interface and phase shift of each 
interface does not occur. If we regard a as 0(1) with respect to //, however, Eq. (|75j) can 
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not satisfy the boundary condition Eq. (JBTj) . On the other hand, if we substitute Eq. (J73J) 
into Eq. (fTIj) . it is found that the amplitude and the phase of the liquid-air surface depends 
on the wavelength of a fluctuation of the solid-liquid interface because of the restoring forces 
a. 

By rewriting the second equation of Eq. (|38|) in the dimensional form, 

dib' 

v ' = - Uo h — = -iku f (y)( k expert + ikx], (76) 
and by comparing it with Eq. (fTT)|) we obtain 

Vk = -iku f(y)(k, (77) 

where f(y) is given by Eq. (J73)) in the long wavelength approximation retaining up to the 
first order in fi . 

IV. GENERAL SOLUTION FOR THE PERTURBED TEMPERATURE DISTRI- 
BUTION IN THE LIQUID 

In the preceding section, we have determined U(y) and v k in Eq. (fTTj) . Next, we must 
determine the amplitude of the perturbed temperature in the liquid under this primary 
shear flow and the amplitude of perturbed normal flow. Since the Peclet number Vh /K>i 
associated with the crystal growth velocity is very small when the typical values of V ~ 10~ 6 
m/s, h ~ 10" 4 m in the experiments 0, 10 1, and Ki = 1.3 x 1CT 7 m 2 /s of water are used, 



we can neglect the second term of Eq. © , then the solution is 

fi{y) =T m - G lV , (78) 

where Gi = (T m — TJ a )//i is unperturbed temperature gradient in the liquid, T\ a is the 
temperature of the liquid-air surface. If we make the substitutions of y = h (l — z), [i = kho 
and UqIio/ki = Pe, which is the Peclet number associated with the flow velocity at the free 
surface, into Eq. (fTTj) . we obtain 

d 2 9i J 2 / pVho \ 2 oh\ . \ . 2 

^ i ^ 1 J ir 1 i 91 91 

= ^Pef(z)expS.-^(l-z)\G l C k , (79) 
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where we have used Eqs. (JUSI, (|77j) and (|75|) . When we put the right hand side of Eq. (J7§J) 
equal to zero, Eq. (J79)l becomes to the equation for a parabolic cylinder function, 

d 2 (f) J 2 / pVh \ 2 ah 



dz 2 



~ Y + (^r) + + ^ Pe 1 + ^ Pe/0 = °" (80) 



Using the confluent hypergeometric function 1F1, the general solutions of Eq. (jSUJ) are 
given by [3| 



ex P (-iHMP,)^) ^ g {l + , i, (^Pe)Vv) , (81) 



= ^exp ^--H>Pe)VVj 1 F 1 ^- + - |1 + — p - 1/2 |, -, (-vtfer'VJ . (82) 
Then we can show that the Wronskian W of the two solutions <pi(z) and 02 (2;) becomes 

W(z)=Mz)^-M*) d -^ = l- (83) 
Therefore, the solution of Eq. (|75|) is given as follows: 

0(2) = + B 2 Mz) + ^Pe / ~ M*)<hV)} fiz'WGta, (84) 



where B\ and -B2 are constants with respect to z, and in Eq. (J84|) we have omitted the 
exponential term on the right hand side of Eq. (f79|) because Vh /ni < 1. In Eqs. (|5T|) 
and (|82|1. we have omitted the terms (pVho/2Ki) 2 and ah^/ni in Eq. ()80|) because we can 
evaluate the ratio of the second term to the first one, (pVh /2Ki) 2 / p 2 = (pV/2nik) 2 <C 1, 
and the ratio of the third term to the first one, <t/iq/(k;/x 2 ) = a/{nik 2 ) <C 1. We are 
concerned with the wave number region that satisfy the former condition. While the latter 
condition amounts to neglecting the time dependence of the perturbed temperature field. 
It was assumed from the outset in the O-F model. This will be justified in Sec. VI. The 
constants B\ and B 2 must be determined from the boundary conditions at the liquid-air 
surface. 

The equation for the temperature distribution in the air is 

v oj^^j**) M (85) 



dt dy \dx 2 dy 
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where T a and T' a are unperturbed temperature and perturbed temperature of the air re- 
spectively, and K a is the thermal diffusivity of the air. The solution for the unperturbed 
temperature field is 

y ~ ho 



T a (y) = ^ + (Ti a - T^) exp 



L 



(86) 



where Tx, is the temperature of the air at y = oo and l a = n a jV is the thermal diffusion 
length of the air. Suppose that the perturbed temperature distribution of the air is expressed 
in the following form: 

V 



T 'a = 9a(y)exp 



2n a 



(y - h ) 



exp[crt + ikx], 



(87) 



where 



9a(y) = T ka exp[-q(y - ho)], (88) 

and T ka is the amplitude of the perturbed temperature of the air. Substitution of Eq. (JHTj) 
into Eq. (|85j) gives 

q = Jtf+ ( +—. (89) 



2,K a I K n 



In the quasi-stationary approximation a/(K a k 2 ) <C 1 and kl a ^> 1, we can approximate 
q = k. 

The continuity of the temperature at the liquid-air surface, y — x) = h + exp [at + 
ikx], is 

(f l + TX = s = (f a + TX=t = T la . (90) 
Linearizing Eq. (|50|l at y — h Q , Eq. (j9T!|) gives to the zeroth order in £ fc , 



1 I 7^ I 

i | y=/io <L\y=ho A lai 



(91) 



and to the first order in 



-Giik + 9i\ y =h exp 



' 2«, 



-C a ^ + T fca — 0, 



(92) 



where G a = (T a — Too)// a - Hereafter we omit the terms including Vho/ni because Vho/ni <C 
1. Heat conservation at the liquid-air surface is 



-K, 



djft + T{) 
dy 



-K, 



d(f a + T' a ) 
dy 



(93) 
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where K a is the thermal conductivity of the air. Similarly, linearizing Eq. (J9H)) at y = ho, 
Eq. (|9*3)) gives to the zeroth order in 

tf,G, = K a G a , (94) 

and to the first order in 

tf,B 2 = /iK a T fca . (95) 
From the first equation of Eq. (|92|L we obtain 

Bi = G,& = -J|,= G z a, (96) 

where we have used the relation Eq. (J71j) in the dimensional form. Here f(z) has the 
following form by substitution of = 1 — z into Eq. (|75j): 

/(z) = (— 6 + zaz + 6z 2 — zaz 3 ) 

6 — 2« 

-2l0&i 144 + <- 174 + 5!a ) 2 - 144 " 

+(210 - llia)z 3 + (-42 + 7ia)z 5 + (6 - za)z 7 } . (97) 

From the second equation of Eq. (|92j). 

Tfca = G a £fc. (98) 
Eliminating T fca from Eqs. (j93|) and (J35j) and using Eq. (|94p. we obtain 

B 2 = A^i- (99) 
By substitutions of Eqs. ()96j) and ()99|) into Eq. (|84|). we finally obtain 



/U=o + VH*)) + ^Pe / {<h(*)MJ) - 0i(^)0 2 (^')} /(^' 



G^Cfc 

#i(*)G,Ck. (100) 



V. APPLICATION 

In this section, we apply the solutions obtained in Sec. Ill and IV to the general formulae 
in Sec. II. If we assume T s = T m in the solid, Eq. (|T9"j) becomes 

LV = KiGi. (101) 
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We solve Eq. (|12|) in the quasi-stationary approximation a/(K. 9 k 2 ) <C 1 and kl s ^> 1, where 
l s = k s /V being the thermal diffusion length of the solid and in the condition that the 
disturbance must vanish far from the solid-liquid interface, then the propagator in the solid 
is 



Qs 



dg s 
dy 



y=o 



9s\y=0 



k. 



(102) 



If we are interested in the long wavelength region such that Qi^i/V, Q s l s ^> 1, using Eq. 
(PUD , Eq. (|2I|> becomes 

G(k) 



, = V Ql (l + ^) + nV Qs ^ 



T? dHl 



where 



and 



and where 



Hi\ z=1 = -/|*= o (0i|*=i + H4>2\z=i) + ifiPeI\ z=1 , 



dH, 



dz 



2=1 



~f\z=0 



dz 



z=i dz 



z=l 



+ i/j,PeJ\ z= i, 



Hz) 



(103) 

(104) 
(105) 
(106) 



which describes the disturbance of the steady state temperature distribution by fluid flow 
normal to the interface, and 



J(z) 



dz 



dz 



(107) 



In the absence of flow, we put Pe = in Eqs. (JHTj) and (f82"|). If we expand Eqs. (jSTjl and 
(182)1 with respect to the powers of \x up to infinity using the recursion relation by setting 
ai = /x 2 /2 and 02 = (see Appendix A), we obtain at z — 1 



>l\z=l 



[i 2 /1 4 



(108) 



>2\z=l 



n fi 2 /i 4 



In the same way, the derivative of Eqs. (J8T|) and (|52*|) at z — 1 are 

d<Pi ( U 3 U 5 

^^ 1 = /X l /X+ ¥ + 5! + - 



(109) 



(110) 
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Then the propagator in the liquid becomes 



-1 dHi 
L dz 


d<h 

z=l 1 dz 


2 = 1 


2 = 1 


h H t 


2=1 ho 4>i 


2 = 1 + /i02 2 = 1 



Qi 



If we take G(k) in Eq. ()103|) as the Gibbs-Thomson effect 

G(k) = -^-k\ 



k. 



112) 



0: 



113) 



where V is the solid-liquid interface tension, Eq. 
the Mullins-Sekerka theory 



reduces to the dispersion relation in 



o T = Vk 



1 - d ^(l + n)k 2 



(7; = 0, 



;ii4) 



where n = K s /Ki and do = T m TC p / L 2 is the capillary length, C p being the specific heat at 
constant pressure, and from Eq. (|l()lj) . 

G p Kl(T m — Ti a ) 



V 



Lhn 



(115) 



It should be noted that the k term in front of the brackets in Eq. (|114|) comes from Eqs. 
(Il()2jl and ()112j) . The \i term in Eq. ()112|) appears as the result of the assumption that the 
heat in the air is transported by thermal diffusion. 

On the other hand, in the presence of flow, in the long wavelength region of about 1cm 
which is the typical wavelength of the wavy pattern observed on the surface of icicles, we 
can neglect the Gibbs-Thomson effect as discussed in the article [8j. From Eq. |T6|) . 

0,U= ( 1 + ^)^C fc , (116) 
and noting gi\ z =\ = Hi\ z= xGiCk from Eq. (jlOOj) . the following relation must be satisfied: 

G(k) = (H l \ z=1 -l)G l . (117) 



Then Eq. (jl(J3j) can be written as 

V [dHi 



a 



h n dz 



2 = 1 



+ Tift (Hi\ z=1 - 1] 



(118) 
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We note that the second term on the right hand side of Eq. ()118|) represents the thermal 
diffusion of latent heat produced by a disturbed solid-liquid interface into the solid. In the 
long wavelength region, we can make approximation of neglecting the fi 2 term in Eqs. (|81|) 
and ()82|). This term is originated from the diffusion term d 2 T[/dx 2 in Eq. (jSJ). In this 
case, the transport of heat in the liquid is dominated by shear flow. Noting that /iPe = 
uohyk/ ' Ki ~ 0(1) for the wavelength of about 1 cm observed on the surface of icicles when the 



typical values of uq ~ 10 2 m/s, ho ~ 10 4 m in the experiments 



lfl and Ki = 1.3 x 10" 7 



m 2 /s of water are used jq, we expand Eqs. (|HT|) and (JHU) with respect to the powers of 
fjPe up to the second order using the recursion relation by setting a\ = (/iPe) 1//2 / (2^/2) and 
a 2 = /iPe/2 - ( / uPe) 1/2 /(2 v / 2) (see Appendix A) as follows: 



1+ -— z A + — ; 
V 24 360 



— z s ) (/xPe) 2 + i(-z 2 - -z 4 ) /iPe, 
Q72 y \2 12 ' ' 



(119) 



z + 



1 



120 2520 

We evaluate each function and its derivative at z — 1 



1 + z — uPe 
12 P 



239 
10080 



(AiPe) 



(120) 



121} 



7 13 
y- = 1 + i 60" Pe - 3360(" Pe ) 2 ' 



122) 



2 13 



(123) 



f/2 



^-l+^Pe-^Pe) 2 , 



-6 24/iRea 



(124) 



(125) 



and we evaluate Eq. (|9"Tj) at z = 0: 

^^ _0 6 — ia 35(6 — 2a) 2 
Substitutions of Eq. (jHZD and Eqs. (ITTD1 to (TT241 into Eqs. (ITUoT) and (ITU7I) and integra- 
tion of them give respectively at z = 1, 



2/iPeJ 



2 = 1 



36 + a 2 
+i 



-a(uPe) + ( — + -^-a 2 ) (uPe) 2 
5 VP 7 V280 3360 J VP ; 

15 + ^a 2 ) (/iPe) + ^aGuPe) 2 



(126) 
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i//PeJ| 



1 



2 = 1 



a(uPe) + I h 

2 VP ; V 35 1440 



7S 
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a" 



36 + a 2 

i i - ( 24 + ^fa a ) (//Pe) + ^a(/iPe) 



(/iPe) 



101 



336 



(127) 



where we have carried out integration by neglect of the first order term in fi in Eq. ()97|) 
because this term is expected to give very small correction to Eqs. ()126|) and ()127|) . By 
substitutions of Eqs. ()121|) to ()127|) into Eq. ()118j) . the final forms of a r and v p = —Oijk for 
the fluctuation of the solid-liquid interface in the long wavelength region are 



Or 



V_ 

h () 



-§a(//Pe) + fj, {36 - §a(/iPe) } -^a(//Pe) + \i {36 - ^a(/iPe) } - 



36 + a 2 



36 + a 2 



28) 



Vr, = 



-±a 2 (/xPe) + /x {6a + 9(/xPe)} 
36 + a 2 



n/i 



6a - ^a 2 (/xPe) + // {6a + f (/xPe)} 
36 + a 2 



(129) 

where we have neglected the first order term in /i in Eq. ()125|) by the same reason as 
mentioned above. Although we have expanded <p\ and 02 up to the second order with 
respect to /zPe as in Eqs. (|119|) and (|120j) . the values of coefficients of (/iPe) 2 are very small 
compared to those of /zPe. Indeed, we have confirmed that the form of a r and v p including 
up to (/iPe) 2 is almost the same as Eqs. ()128j) and ()129j) in the long wavelength region such 
that k < 10 3 / m. Therefore, it is sufficient to approximate a r and v p up to the first order in 
/iPe. 

The rate of volume flow down the inclined plane in the experiment [l0| is 



Q = u Q l 



ho 



n v y i , 



'-u h l, 



(130) 



where I is the width of the gutter, and 



9 h o ■ a 
u = sin 6 

2v 



;i3i) 



is the surface velocity Q|. If we eliminate Uq from Eqs. (|130|) and (|131|) . the mean thickness 
ho of the liquid can be expressed with respect to Q and 6: 



3uQ 
gl sin 9 



1/3 



(132) 
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Then, pPe and a can be expressed in terms of ho, respectively, 

gsm9 4 



/iPe = »- h\k, (133) 



a = 2 cot 9h k + a 2 h k 3 , (134) 
where we have defined the capillary constant associated with the surface tension 7 of the 

n 

liquid-air [13j : 

(135) 



gpi sin 9 

We note that this capillary constant depends on 9, and that this typical value is about 3.9 
mm for 7 = 7.6 x 10 -2 N/m and pi = 1.0 x 10 3 Kg/m 3 of water when 9 = 7r/2. 

From Fig. El to H we use the values of T la = -0.06° C, Q = 160 ml/hr and I = 0.03 
m in the experiments , and the physical properties of water, L = 3.3 x 10 8 J/m 3 , 

C p = 4.2 x 10 6 J/(Km 3 ), «, = 1.3 x 10" 7 m 2 /s, v = 1.8 x 10~ 6 m 2 /s, 7 = 7.6 x 10~ 2 N/m 
and n = K s /Ki = 3.92, where K s is the thermal conductivity of ice. The reason for choosing 
the value of Q = 160 ml/hr is that the clearest wavy pattern was observed at this value in 
the experiment [10 . Since the crystal growth velocity V observed in the actual experiment 
is about 10 -6 m/s 9j, from Eq. (|115|) we obtain the value of T\ a = —0.06° C for water when 
ho = 10 _4 m. Although T} a is to be determined by the condition of the surrounding air, we 
use this value for T\ a to determine the value of V from Eq. pi 5)1 when varing 9. 

The solid line in Fig. |21 shows the amplification rate Eq. (|128|) versus wave number k 
for 9 = 7r/2. This shows that a r takes a maximum value (T rmax at a wave number. The 
characteristic time for most unstable mode is l/cr rmax and is about 30 minutes in this case. 
Indeed, it is reported in the experiment that a periodic structure as the original form of 
wavy pattern is observed in about 30 minutes |l0j. On the other hand, the dashed line in 
Fig. El shows a r when we neglect the restoring forces due to gravity and surface tension, i.e., 
when we put a = in Eq. ()128|) . Then a r is always positive in the range of our interests. 

Fig. El shows the phase velocity Eq. ()129|) versus wave number k for 9 = tt/2. This shows 
that the fluctuation of the solid-liquid interface for the maximum point of a r in Fig. El m oves 
upward along the icicle with the magnitude of about 0.6 V. Indeed, there is evidence to 
support our predictions that many tiny air bubbles dissolved in the thin flowing liquid are 
trapped in just upstream region of any protrudent part on a growing icicle and its region 
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a 




FIG. 2: The amplification rate oy versus wave number k for T/ a = —0.06° C, Q = 160 ml/hr and 
6 = 7r/2. Solid line : with restoring forces. Dashed line : without restoring forces. 




FIG. 3: The phase velocity v p = —<Ji/k versus wave number k for T} a = —0.06° C, Q = 160 ml/hr 
and 9 = vr/2. 

n 

migrates in the upward direction during growth (see Fig. 9B in Ref. [9]). This suggests that 
the velocity of ice growth is faster in the upstream region of each protrudent part. On the 
other hand, in the O-F model, it was predicted that the fluctuation moves downward along 
the icicle with the phase velocity, 

v p = V ^— = 136 

{l-lilo(^) 2 } 2 + {^Pe} 2 

If this prediction is correct, air bubbles would be trapped in the downstream region of each 
protrudent part and migrate in the downward direction during growth. 

Fig. |H shows the dependence of the wavelength \ max obtained theoretically or the mean 
wavelength X mean obtained in the experiment on the angle 9 of the inclined plane. The 
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closed squares represent mean wavelength obtained by the experiment 



0.83 

mean ~ (sinfl)°-6~°-9 (Cm) ' ( } 

Using Eq. (|128j) . we determine the wavelength A max at the maximum point of ay, and its 
dependence on the angle is found to be 

0.98 , , , , 

Xmax ~ (sinfl)0-6~0-65 ( cm )' ( 138 ) 

which is shown in the closed circles. We note that this dependence of A max on sin# comes 
from not only /iPe of Eq. f!133|) but also cot#, h and a in Eq. (j!34|) . Our results are in 
good agreement with experiment. On the other hand, the amplification rate obtained in the 
O-F model is 
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(//Pe) 



a r = Vk 10080 V ' (139) 

{l-lilo(^) 2 } 2 +{^Pe} 2 

The closed triangles are X max at the maximum point of Eq. (|139|) . Then the result is 

0.47 

Xmax ~ (sing)V3 ( cm )' ( 14 °) 
In this case, we note that the dependence of X max on sin (9 comes from only fiPe of Eq. f!133|) . 



VI. DISCUSSION 

Some differences between our results and their results jsj appear to arise from the following 
reasons. A main difference is originated from the order estimate of Eq. (jUHj) or Eq. ()134|) . 
If we use the values of h Q at Q = 160 ml/hr and k = 27r/\ max , where X ma x is taken from 
Eq. ()138j) . the values of a and fi for 0.1 < sin# < 1 take a range, 0.4 < a < 0.8 and 
0.03 < < 0.06, respectively. Therefore, we have treated a as zeroth order in terms of 
fi. On the other hand, it was regarded as first order in [i in the O-F model. When we 
determine the perturbed stream function, these differences cause different forms between 
Eqs. (|73^1 and (f?3j) . As a result of that, different dependence of A max on sin 6 between Eqs. 
(|138|) and (|140j) has occured. The sin# term in Eq. (|128|) is included not only in fiPe but 
also a. On the other hand, in the O-F model, the sin# term in Eq. ()139j) appears in only 
/xPe. If we use Eq. (J1HJ), in which the effect of restoring forces due to gravity and surface 
tension on the liquid-air surface is included, a r takes a maximum value at a wave number. 
On the other hand, if we use Eq. (J7SJ), in which the effect of restoring forces is not included, 
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FIG. 4: The dependence of X ma x or A mean on the angle of inclined plane for Q = 160 ml/hr. 



Closed circles : present result. Closed triangles : Ogawa-Furukawa's result 
experimental result 



. Closed squares : 



a r is always positive in the long wavelength region. In spite of absence of the a term in 
;he O-F model, the similar curve as solid line in Fig. El was obtained (see Fig. 4 in Ref . 
^I). The existence of maximum of their a r is the result of expansion of the temperature 
fluctuation in the liquid up to (/iPe) 2 , for example, which is reflected in the numerator in 
Eq. (j!39|) . However, we have confirmed that it is sufficient to approximate a T up to the first 
order in /zPe. The existence of maximum of our a r comes from the effect of a and /zPe. 

These different results and the different prediction of the direction of the phase velocity 
between Eqs. ()129|) and (jl36|) may be due to the difference of the boundary condition for 
the temperature at the solid-liquid interface and that of the liquid-air surface between ours 
and theirs. Instead of Eqs. (fE?j) and (J§Uj) . the following boundary conditions were used in 
the O-F model, respectively, 

(T, + Tl)\ y=( = (f s + T%=C = T m , (141) 

(T l+ T{)\ y= s = (f a + T>)\ y =t- (142) 
In order to determine the two constants B\ and B2 in Eq. (|H4"|) independently, and in the 
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absence of flow, to recover the usual Mullins-Sekerka theory from the general dispersion 
relation Eq. ([103)1 . we have used the boundary condition Eq. (|9Uj) instead of Eq. ()142|) . For 
long wavelength fluctuation of about 1 cm of the solid-liquid interface, since we can neglect 
the change of the melting temperature due to the Gibbs-Thomson effect, adopting Eq. (|141|) 
appears to be appropriate. In the presence of flow, however, Eq. ([11 7|) suggests that there 
exists a shift of the melting temperature depending on the wave number. Therefore we have 
used the boundary condition Eq. (fTB"|) instead of Eq. ()141[) . 

/x term in the numerator in Eq. ()128|) is the cause of instability. This term is originated 
from spatial derivative of the perturbed air temperature distribution at the deformed liquid- 
air surface as indicated on the right hand side of Eq. ()95|) . From Eq. (|68|) or Eq. (|134|) , since 
the value of a is very small in the long wavelength region, the effect of restoring forces due 
to gravity and surface tension on the liquid- air surface is small. Therefore, in the low wave 
number region instability effect with positive terms in Eq. (|128|) dominates stability effect 
with negative terms. On the other hand, as increasing the wave number, since the value of 
a increases, the effect of restoring forces is large. Then a(/iPe) and a 2 terms with negative 
sign in the numerator of Eq. ()128j) dominate the instability terms. As a result of that, we 
obtain the solid curve in Fig. |21 In order to interpret correctly the physical mechanism of 
instability and stability of the solid-liquid interface and why the solid-liquid interface moves 
in the upstream direction, it is necessary to understand the relative phase of modes at each 
interface using the relation Eq. (fT4"j) and a shift of melting temperature due to flow and 
restoring forces as suggested in Eq. ()117j) . This will be shortly clarified in another paper. 

In Sec. Ill and IV, we have made some assumptions. Here we justify them. We have 
assumed the time independence of the purterbed flow, therefore we have neglected the a*/ ' ji 
term in Eqs. (jH3), (p^j) . (p>4"|) and (jo^j) . This assumption is valid because we see from Fig. |2 
and Fig. |3]that ~ 10~ 6 for V ~ 10 -6 m/s. Therefore the condition cr*/// <C 1 is satisfied. 
The same can be said for the fluctuation of the temperature. For a deformation of wave 
number k, the characteristic delay time of the fluctuation of the temperature is of order 
^thermal ~ (/fyA; 2 ) -1 . /St t hermai is much smaller than the characteristic time of evolution of 
the mode, a~ l . Therefore the condition a r /(K,ik 2 ) <C 1 is satisfied. These mean that the 
perturbed flow field and the perturbed temperature field respond relatively rapidly to the 
slow development of the solid-liquid interface. 
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VII. CONCLUSION 



The restoring forces due to gravity and surface tension determine the shape of free surface 
and do not directly act on the solid-liquid interface. However, the effect of restoring forces 
has played an important role on stabilization of the solid-liquid interface. Although the 
Gibbs-Thomson effect acts effectivley on the micrometer scale, we have found that the effect 
of restoring forces is more effective for long wavelength fluctuation of the order of mm, which 
is of the order of the capillary constant associated with the surface tension of the liquid-air. 
Therefore the wavy pattern observed on the surface of icicles and inclined plane occures 
on longer length scales compared to the wavelength predicted by the usual Mullins-Sekerka 
theory. 

Since our calculations are based on the linear stability analysis, our formulations and 
the Ogawa and Furukawa's formulations do not have direct correspondences. The relation 
between our present formulations and previous ones is under investigation by Ogawa and 
Furukawa. 
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APPENDIX A: 

We separate Eq. (J81|) into the real part and imaginary part: 
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Using the first two coefficients, 



ai = ^ 2 + 27! (/uPe)1/2 ' 



(A2) 



a2 = ^ Pe "272 (/iPe)1/2 ' 



(A3) 



the other coefficients a, for odd numbers are obtained from the following recursion relation, 



tt2j+1 = (j + l)(2j + 1) l® 1 ^' 1 ~ a2a2j + ^(^ Pe ) 1/2 ( a 2i-i + a 2i) } U = 1, 2, 3, . . .) 



(A4) 



and for even numbers, 
1 



0-2j+2 



(j + 1)(2j + 1) 



2j 



+ a 2 a 2j -i -=(/iPe) 1/2 (a 2 j-i - o^-) > (j = 1, 2, 3, . . .). 



v/2 



(A5) 



Next, we separate the derivative of Eq. (jHTjl into the real part and imaginary part: 

d(j>x{z) 



Using 



(/iPe) 1/2 2 
(/iPe) 1 / 2 



cos 



— sin 



2^ 

, (^Pe) 1/2 . 



j=0 



(EW* +1 : 



2^ 



j=0 



61 = 2ai - ^(AiPe) 1/2 , 
v2 



(A6) 



(A7) 



6 2 = 2a 2 + — (^Pe) 1/2 , 
v2 



the other coefficients 6, for odd numbers are obtained from 



hj+i = 2(j + l)a 2i+ i 
and for even numbers, 

b 2 j+2 = 2(j + l)a 2j+2 



V2 



(jjPe) 1 P{a 2j - 1 + 023) {] = 1,2,3,...), 



-J=(/iPe) 1/2 (a 2i _ 1 -a 2j ) 



(j = 1,2,3, 



(A8) 



(A9) 



(A10) 
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Eq. (|82|) is separated into the real part and imaginary part in the same way: 



Mz) = exp 



Using 



j („Pe)'/ 2 A [ / ( P Pe)V 2 V f 2j+3 



sin 



(/xPe) 1 / 2 
2\/2 



)oo 
i=o 



(/iPe) 1 / 2 ,\ ^ 



uPe W 2 , 
+ sin | VP i. 2- 
2^ 



j=0 

OO 



' 3=0 



-1^' +3 ) 



(All) 



(A12) 



c 2 = \a 2 - ^(^e) 1/2 , 
the other coefficients Cj for odd numbers are obtained from 

i r 2 j + 1 



C2j+1 



0' + l)(2j + 3) 



a l c 2i'-l — fl2 c 2i + 



(/iPe) 1/2 (c 2j _i + c 2i ; 



and for even numbers, 
1 

C2J+2 " (j + l)(2 J + 3) 



aic 2 j + a 2 c 2 



(A13) 



(j = 1,2,3,...), 
(A14) 



i-i ~ ^^(/iPe) 1 / 2 (c 2j _ 1 - c 2i )| (j = 1, 2, 3, . . .). 

(A15) 



Finally, derivative of 02 is separated into the real part and imaginary part: 



(/xPe) 1 / 2 2 

exp = — z 

■ 1 2^ 



cos 



— sin 



(/xPe) 1 / 2 
2v"2 



(EW +a ) 

j=0 



+.<COS(^^, 2 )(^^ + 2^ +2 ) 

3=0 



+ sin 



2^2 



2j+2^ 



3=0 



where c?i = a 1; d 2 = a 2 , and other coefficients eL for odd numbers are obtained from 



(2j + 3)c 2i+1 - -J=(/iPe) 1 / 2 (c 2j _i + c 2j ) (j = 1, 2, 3, 



(A16) 



(A17) 
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and for even numbers, 

d 2j+2 = (2j + 3)c 2i+2 - -^=(ftPe) 1/2 (c 2j ^ - c 2i ) (j = 1, 2, 3, . . .)• (A18) 

It should be noted that all coefficients are obtained from only a% and a 2 . Eqs. (j!21|) to 
()124|) are valid only in the long wavelength region because we neglect the /i 2 term in Eq. 
(IA2J) . This means that heat transport is dominated by shear flow. On the other hand, in 
the absence of flow, we put Pe = in Eq. ()A2jl . Then Eqs. (jl()8|) to are obtained. In 

this case, heat transport is dominated by thermal diffusion. 
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